altKrigCRF<-function (krig.x, krig.y, resid.x, resid.y, resid.direction, 
    Model, Nugget = 0, Range, sill, Smooth = FALSE, bandwidth, 
    Plot = FALSE, Xlim = NULL, Ylim = NULL, PlotVar = FALSE, 
    ...) 
{
    if ((length(krig.x) != length(krig.y)) | (length(resid.x) != 
        length(resid.y)) | (length(resid.x) != length(resid.direction)) | 
        (length(resid.y) != length(resid.direction))) 
        stop("lengths of vector inputs unequal")
    if ((Nugget < 0) | (Nugget > 1)) 
        stop("Nugget invalid")
    xx <- sort(krig.x)
    yy <- sort(krig.y)
    nx <- length(xx)
    ny <- length(yy)
    krig.y <- yy
    krig.x <- xx
    Distances <- as.matrix(dist(cbind(resid.x, resid.y)))
    Ncol <- ncol(Distances)
    K <- c()
    for (i in 1:Ncol) {
        K <- cbind(K, sill + (1 - Nugget - sill) * CovarianceFct(x = Distances[, 
            i]/Range, model = Model, param = c(mean = 0, variance = 1, 
            nugget = 0, scale = 1, ...), dim = 1, fctcall = "Cov"))
    }
    diag(K) <- 1
    Kinv <- solve(K)
    U <- t(cbind(cos(resid.direction), sin(resid.direction)))
    n <- length(krig.x)
    krig.direction <- vector(mode = "numeric", length = n)
    krig.variance <- krig.direction
    for (i in 1:n) {
        distances <- sqrt((krig.x[i] - resid.x)^2 + (krig.y[i] - 
            resid.y)^2)
        c <- sill + (1 - Nugget - sill) * CovarianceFct(x = distances/Range, 
            model = Model, param = c(mean = 0, variance = 1, 
                nugget = 0, scale = 1, ...), dim = 1, fctcall = "Cov")
        c[distances == 0] <- 1
        w <- Kinv %*% c
        u <- U %*% w
        krig.direction[i] <- atan2(u[2], u[1])
        krig.variance[i] <- 2 - 2 * sqrt(as.numeric(t(c) %*% 
            Kinv %*% c))
    }
    if (Smooth) {
        xx.dx <- xx[2] - xx[1]
        yy.dy <- yy[2] - yy[1]
        ImageList.x <- as.image(cos(krig.direction), x = data.frame(krig.x, 
            krig.y), nrow = nx, ncol = ny, boundary.grid = FALSE)
        smooth.x <- image.smooth(ImageList.x, theta = bandwidth)
        ImageList.y <- as.image(sin(krig.direction), x = data.frame(krig.x, 
            krig.y), nrow = nx, ncol = ny, boundary.grid = FALSE)
        smooth.y <- image.smooth(ImageList.y, theta = bandwidth)
        krig.direction <- as.vector(t(atan2(smooth.y$z, smooth.x$z)))
    }
    if (Plot) {
        if (!PlotVar) {
            plot(krig.x, krig.y, ty = "n", xlab = "", ylab = "", 
                asp = 1, xlim = Xlim, ylim = Ylim)
            arrow.plot(a1 = krig.x, a2 = krig.y, u = cos(krig.direction), 
                v = sin(krig.direction), arrow.ex = 0.06, xpd = FALSE, 
                true.angle = TRUE, length = 0.05, col = 1)
        }
        else {
            krig.variance.matrix <- matrix(data = krig.variance, 
                nrow = nx, ncol = ny, byrow = TRUE)
            if (!is.null(Xlim)) {
                filled.contour(x = xx, y = yy, z = krig.variance.matrix, 
                  xlim = Xlim, ylim = Ylim, color = terrain.colors, 
                  key.title = title(main = "Circ Krig \nVariance", 
                    cex.main = 0.8), asp = 1, plot.axes = {
                    axis(1)
                    axis(2)
                    points(resid.x, resid.y, pch = 20, cex = 0.65)
                  })
            }
            else {
                filled.contour(x = xx, y = yy, z = krig.variance.matrix, 
                  color = terrain.colors, key.title = title(main = "Circ Krig \nVariance", 
                    cex.main = 0.8), asp = 1, plot.axes = {
                    axis(1)
                    axis(2)
                    points(resid.x, resid.y, pch = 20, cex = 0.65)
                  })
            }
        }
    }
    else return(list(x = krig.x, y = krig.y, direction = krig.direction, 
        variance = krig.variance))
}
